#ifndef __CXFLOAT__
#define __CXFLOAT__

#include "refloat.h"
#include <vector>
#include <algorithm>

class cxfloat
{
public:
  cxfloat();
  cxfloat(std::nullptr_t ptr);
  cxfloat(long int re, long int im, mp_prec_t prec=4, bool temporary=false);
  cxfloat(std::string re, std::string im, mp_prec_t prec=4, int base=16, bool temporary=false);
  cxfloat(const cxfloat& z);
  void clear();
  void clear_real();
  void clear_imag();
  ~cxfloat();
  void set_temporary(mpfr_prec_t prec);
  void set_temporary_real(mpfr_prec_t prec);
  void set_temporary_imag(mpfr_prec_t prec);
  void set(const cxfloat& x);
  void set(const refloat& x, const refloat& y);
  void set(long int re, long int im, mp_prec_t prec);
  void set(std::string re, std::string im, mp_prec_t prec, int base);
  void set(long int re, long int im);
  void set(std::string re, std::string im, int base);
  void set(std::nullptr_t ptr);
  void point( const cxfloat& z );
  void set_imag(long int im);
  void set_imag(const cxfloat& x);
  void set_real(long int im);
  void set_real(const cxfloat& x);
  void set_prec(mpfr_prec_t prec);
  mpfr_prec_t get_prec() const;
  void neg();
  void add(const cxfloat& x);
  void add(const refloat& x);
  void sub(const cxfloat& x);
  void sub(const refloat& x);
  void mul(const cxfloat& x);
  void mul(const refloat& x);
  void sqr();
  void div(const cxfloat& x);
  void div(const refloat& x);
  void sqrt();
  void mul_power_2(long int p);
  void mul_by_i();
  void conj();
  void reciprocal();
  void print(std::string &re, std::string &im, int base=16) const;
  void print_real(std::string &re, int base=16, mpfr_prec_t prec=0, bool bold=true) const;
  void print_imag(std::string &im, int base=16, mpfr_prec_t prec=0, bool bold=true) const;
  void print_mantissa(std::string &re, std::string &im, int base=16) const;
  void get_real(refloat& val) const;
  void get_imag(refloat& val) const;
  bool isTemp() const;
  void prefetch() const;
  
  cxfloat& operator =(const cxfloat& b);
  
  __mpfr_struct* real_part;
  __mpfr_struct* imag_part;
  
private:
  bool is_temp;
};

typedef std::vector<cxfloat> cxfloat_vec;


extern thread_local cxfloat *cxfloat_cache1;
extern thread_local cxfloat *cxfloat_cache2;
extern thread_local bool cxfloat_all_allocations_temporary;
extern thread_local std::vector<__mpfr_struct*>* cxfloat_temporaries;

inline void setCxAllocationsToTemporary(bool tt){cxfloat_all_allocations_temporary=tt;}


void init_cxfloat_cache(mpfr_prec_t p);
void set_cxfloat_precision(mpfr_prec_t p);
void clear_cxfloat_cache();

void init_cxfloat_cache(void *arg);
void set_cxfloat_precision(void *arg);
void clear_cxfloat_cache(void *arg);




void swap(cxfloat& a, cxfloat& b);
void neg(const cxfloat& a, cxfloat& c);
void add(const cxfloat& a, const cxfloat& b, cxfloat& c);
void sub(const cxfloat& a, const cxfloat& b, cxfloat& c);
void mul(const cxfloat& a, const cxfloat& b, cxfloat& c);
void mul(const refloat& a, const cxfloat& b, cxfloat& c);
void sqr(const cxfloat& a, cxfloat& c);
void div(const cxfloat& a, const cxfloat& b, cxfloat& c);
void sqrt(const cxfloat& a, cxfloat& b);
void mul_power_2(const cxfloat& a, long int p, cxfloat& c);
void mul_by_i(const cxfloat& a, cxfloat& b);
void conj(const cxfloat& a, cxfloat& b);
void magnitude(const cxfloat& a, cxfloat& b);
void magnitude(const cxfloat& a, refloat& b);
void compare_mantissas(const cxfloat& x1, const cxfloat& x2, unsigned long int &d_re, unsigned long int &d_im, int base=16);
void sin_cos(const cxfloat& z, cxfloat& sn, cxfloat& cs);
void cos(const cxfloat& z, cxfloat& cs);
void sin(const cxfloat& z, cxfloat& sn);
void exp(const cxfloat& z, cxfloat& e);
void log(const cxfloat& z, cxfloat& l);
void pow(unsigned long int n, const cxfloat& x, cxfloat& v);
int compare_cxfloat_real(const void *a, const void *b);
void sort_by_real(cxfloat *arr, long int size);
int compare_real(const cxfloat& a, const cxfloat& b);
int compare_real_abs(const cxfloat& a, const cxfloat& b);
int compare_imag(const cxfloat& a, const cxfloat& b);
int compare_imag_abs(const cxfloat& a, const cxfloat& b);


inline void cxfloat::prefetch() const
{
  unsigned long int word_size = 8*sizeof(unsigned long int);
  unsigned long int nwords = mpfr_get_prec(real_part)/word_size;
  for(unsigned long int i=0;i<nwords;++i)
  {
    __builtin_prefetch(&((real_part->_mpfr_d)[i]));
  }
  nwords = mpfr_get_prec(imag_part)/word_size;
  for(unsigned long int i=0;i<nwords;++i)
  {
    __builtin_prefetch(&((imag_part->_mpfr_d)[i]));
  }
}


inline cxfloat& cxfloat::operator =(const cxfloat& b)
{
  if( (real_part == nullptr) && (b.real_part != nullptr) )
  {
    allocate_mpfr_struct(real_part);
    init_mpfr_value(real_part, mpfr_get_prec(b.real_part));
    mpfr_set(real_part, b.real_part, GMP_RNDN);
  }
  else if( (b.real_part == nullptr) && (real_part != nullptr) )
  {
    free(real_part->_mpfr_d);free( real_part );real_part = nullptr;
  }
  else if( (b.real_part != nullptr) && (real_part != nullptr) )
  {
    if( mpfr_get_prec(b.real_part) != mpfr_get_prec(real_part) )
    {
      free(real_part->_mpfr_d);free( real_part );
      allocate_mpfr_struct(real_part);
      init_mpfr_value(real_part, mpfr_get_prec(b.real_part));
    }
    mpfr_set(real_part, b.real_part, GMP_RNDN);
  }
  
  if( (imag_part == nullptr) && (b.imag_part != nullptr) )
  {
    allocate_mpfr_struct(imag_part);
    init_mpfr_value(imag_part, mpfr_get_prec(b.imag_part));
    mpfr_set(imag_part, b.imag_part, GMP_RNDN);
  }
  else if( (b.imag_part == nullptr) && (imag_part != nullptr) )
  {
    free(imag_part->_mpfr_d);free( imag_part );imag_part = nullptr;
  }
  else if( (b.imag_part != nullptr) && (imag_part != nullptr) )
  {
    if( mpfr_get_prec(b.imag_part) != mpfr_get_prec(imag_part) )
    {
      free(imag_part->_mpfr_d);free( imag_part );
      allocate_mpfr_struct(imag_part);
      init_mpfr_value(imag_part, mpfr_get_prec(b.imag_part));
    }
    mpfr_set(imag_part, b.imag_part, GMP_RNDN);
  }
  
  return *this;
}


inline cxfloat::cxfloat() : real_part(nullptr), imag_part(nullptr), is_temp(cxfloat_all_allocations_temporary)
{
  mpfr_prec_t prec = mpfr_get_prec(cxfloat_cache1->real_part);
  if(is_temp==false)
  {
    allocate_mpfr_struct(real_part);
    init_mpfr_value(real_part, prec);
    
    allocate_mpfr_struct(imag_part);
    init_mpfr_value(imag_part, prec);
  }
  else
  {
    set_temporary(prec);
  }
}


inline cxfloat::cxfloat(std::nullptr_t ptr) : real_part(nullptr), imag_part(nullptr), is_temp(cxfloat_all_allocations_temporary)
{
  
}


inline cxfloat::cxfloat(long int re, long int im, mp_prec_t prec, bool temporary) : is_temp(temporary || cxfloat_all_allocations_temporary)
{
  if(is_temp==false)
  {
    allocate_mpfr_struct(real_part);
    init_mpfr_value(real_part, prec);
    
    allocate_mpfr_struct(imag_part);
    init_mpfr_value(imag_part, prec);
  }
  else
  {
    set_temporary(prec);
  }
  mpfr_set_si(real_part, re, GMP_RNDN);
  mpfr_set_si(imag_part, im, GMP_RNDN);
}


inline cxfloat::cxfloat(const cxfloat& z) : real_part(nullptr), imag_part(nullptr), is_temp(z.isTemp() || cxfloat_all_allocations_temporary)
{
  (*this) = z;
}


inline cxfloat::cxfloat(std::string re, std::string im, mp_prec_t prec, int base, bool temporary)
{
  if(is_temp==false)
  {
    allocate_mpfr_struct(real_part);
    init_mpfr_value(real_part, prec);
    
    allocate_mpfr_struct(imag_part);
    init_mpfr_value(imag_part, prec);
  }
  else
  {
    set_temporary(prec);
  }
  mpfr_set_str(real_part, re.c_str(), base, GMP_RNDN);
  mpfr_set_str(imag_part, im.c_str(), base, GMP_RNDN);
}


inline void cxfloat::clear_imag()
{
  if(imag_part != nullptr)
  {
    if((is_temp==false) || (!cxfloat_temporaries)){ free(imag_part->_mpfr_d);free( imag_part ); }
    else{ cxfloat_temporaries->push_back(imag_part); }
    imag_part = nullptr;
  }
}


inline void cxfloat::clear_real()
{
  if(real_part != nullptr)
  {
    if((is_temp==false) || (!cxfloat_temporaries)){ free(real_part->_mpfr_d);free( real_part ); }
    else{ cxfloat_temporaries->push_back(real_part); }
    real_part = nullptr;
  }
}


inline void cxfloat::clear()
{
  clear_real();
  clear_imag();
}


inline cxfloat::~cxfloat()
{
  this->clear();
}


inline bool cxfloat::isTemp() const
{
  return is_temp;
}


inline void cxfloat::point(const cxfloat& z)
{
  real_part = z.real_part;
  imag_part = z.imag_part;
  is_temp = z.is_temp;
}


inline void cxfloat::set(std::nullptr_t ptr)
{
  real_part = nullptr;
  imag_part = nullptr;
}


inline void cxfloat::set(std::string re, std::string im, mp_prec_t prec, int base)
{
  if(is_temp==false)
  {
    if((real_part != nullptr))
    {
      free(real_part->_mpfr_d);free( real_part );
      allocate_mpfr_struct(real_part);
      init_mpfr_value(real_part, prec);
    }
    if((imag_part != nullptr))
    {
      free(imag_part->_mpfr_d);free( imag_part );
      allocate_mpfr_struct(imag_part);
      init_mpfr_value(imag_part, prec);
    }
  }
  else
  {
    if(real_part != nullptr){ cxfloat_temporaries->push_back(real_part); }
    if(imag_part != nullptr){ cxfloat_temporaries->push_back(imag_part); }
    set_temporary(prec);
  }
  mpfr_set_str(real_part, re.c_str(), base, GMP_RNDN);
  mpfr_set_str(imag_part, im.c_str(), base, GMP_RNDN);
}


inline void cxfloat::set(std::string re, std::string im, int base)
{
  if( (real_part == nullptr) || (imag_part == nullptr) )
  {
    mpfr_prec_t p = mpfr_get_prec(cxfloat_cache1->real_part);
    set( re, im, p, base );
  }
  else
  {
    mpfr_set_str(real_part, re.c_str(), base, GMP_RNDN);
    mpfr_set_str(imag_part, im.c_str(), base, GMP_RNDN);
  }
}


inline void cxfloat::set(long int re, long int im, mp_prec_t prec)
{
  if(is_temp==false)
  {
    if((real_part != nullptr))
    {
      free(real_part->_mpfr_d);free( real_part );
      allocate_mpfr_struct(real_part);
      init_mpfr_value(real_part, prec);
    }
    if((imag_part != nullptr))
    {
      free(imag_part->_mpfr_d);free( imag_part );
      allocate_mpfr_struct(imag_part);
      init_mpfr_value(imag_part, prec);
    }
  }
  else
  {
    if(real_part != nullptr){ cxfloat_temporaries->push_back(real_part); }
    if(imag_part != nullptr){ cxfloat_temporaries->push_back(imag_part); }
    set_temporary(prec);
  }
  mpfr_set_si(real_part, re, GMP_RNDN);
  mpfr_set_si(imag_part, im, GMP_RNDN);
}


inline void cxfloat::set(long int re, long int im)
{
  if( (real_part == nullptr) || (imag_part == nullptr) )
  {
    mpfr_prec_t p = mpfr_get_prec(cxfloat_cache1->real_part);
    set( re, im, p );
  }
  else
  {
    mpfr_set_si(real_part, re, GMP_RNDN);
    mpfr_set_si(imag_part, im, GMP_RNDN);
  }
}


inline void cxfloat::set(const cxfloat& x)
{
  if(this == (&x)){return;}
  (*this) = x;
}


inline void cxfloat::set(const refloat& x, const refloat& y)
{
  if( (real_part == nullptr) && (x.value != nullptr) )
  {
    allocate_mpfr_struct(real_part);
    init_mpfr_value(real_part, mpfr_get_prec(x.value));
    mpfr_set(real_part, x.value, GMP_RNDN);
  }
  else if( (x.value == nullptr) && (real_part != nullptr) )
  {
    free(real_part->_mpfr_d);free( real_part );real_part = nullptr;
  }
  else if( (x.value != nullptr) && (real_part != nullptr) )
  {
    if( mpfr_get_prec(x.value) != mpfr_get_prec(real_part) )
    {
      free(real_part->_mpfr_d);free( real_part );
      allocate_mpfr_struct(real_part);
      init_mpfr_value(real_part, mpfr_get_prec(x.value));
    }
    mpfr_set(real_part, x.value, GMP_RNDN);
  }
  
  if( (imag_part == nullptr) && (y.value != nullptr) )
  {
    allocate_mpfr_struct(imag_part);
    init_mpfr_value(imag_part, mpfr_get_prec(y.value));
    mpfr_set(imag_part, y.value, GMP_RNDN);
  }
  else if( (y.value == nullptr) && (imag_part != nullptr) )
  {
    free(imag_part->_mpfr_d);free( imag_part );imag_part = nullptr;
  }
  else if( (y.value != nullptr) && (imag_part != nullptr) )
  {
    if( mpfr_get_prec(y.value) != mpfr_get_prec(imag_part) )
    {
      free(imag_part->_mpfr_d);free( imag_part );
      allocate_mpfr_struct(imag_part);
      init_mpfr_value(imag_part, mpfr_get_prec(y.value));
    }
    mpfr_set(imag_part, y.value, GMP_RNDN);
  }
}


inline void cxfloat::set_imag(long int im)
{
  if(imag_part == nullptr)
  {
    mpfr_prec_t p = mpfr_get_prec(cxfloat_cache1->imag_part);
    if(is_temp == false)
    {
      allocate_mpfr_struct(imag_part);
      init_mpfr_value(imag_part, p);
    }
    else
    {
      set_temporary_imag(p);
    }
  }
  mpfr_set_si(imag_part, im, GMP_RNDN);
}


inline void cxfloat::set_imag(const cxfloat& x)
{
  if(x.imag_part != nullptr)
  {
    if(imag_part == nullptr)
    {
      mpfr_prec_t p = mpfr_get_prec(x.imag_part);
      if(is_temp == false)
      {
        allocate_mpfr_struct(imag_part);
        init_mpfr_value(imag_part, p);
      }
      else
      {
        set_temporary_imag(p);
      }
    }
    
    mpfr_set(imag_part, x.imag_part, GMP_RNDN);
  }
  else if(imag_part != nullptr)
  {
    free(imag_part->_mpfr_d);free( imag_part );imag_part = nullptr;
  }
}


inline void cxfloat::set_real(long int re)
{
  if(real_part == nullptr)
  {
    mpfr_prec_t p = mpfr_get_prec(cxfloat_cache1->real_part);
    if(is_temp == false)
    {
      allocate_mpfr_struct(real_part);
      init_mpfr_value(real_part, p);
    }
    else
    {
      set_temporary_real(p);
    }
  }
  mpfr_set_si(real_part, re, GMP_RNDN);
}


inline void cxfloat::set_real(const cxfloat& x)
{
  if(x.real_part != nullptr)
  {
    if(real_part == nullptr)
    {
      mpfr_prec_t p = mpfr_get_prec(x.real_part);
      if(is_temp == false)
      {
        allocate_mpfr_struct(real_part);
        init_mpfr_value(real_part, p);
      }
      else
      {
        set_temporary_real(p);
      }
    }
    
    mpfr_set(real_part, x.real_part, GMP_RNDN);
  }
  else if(real_part != nullptr)
  {
    free(real_part->_mpfr_d);free( real_part );real_part = nullptr;
  }
}


inline void cxfloat::set_prec(mp_prec_t prec)
{
  if(real_part != nullptr)
  {
    mpfr_set((*cxfloat_cache1).real_part, real_part, GMP_RNDN);
    free(real_part->_mpfr_d);free( real_part );
    allocate_mpfr_struct(real_part);
    init_mpfr_value(real_part, prec);
    mpfr_set(real_part, (*cxfloat_cache1).real_part, GMP_RNDN);
  }
  if(imag_part != nullptr)
  {
    mpfr_set((*cxfloat_cache1).imag_part, imag_part, GMP_RNDN);
    free(imag_part->_mpfr_d);free( imag_part );
    allocate_mpfr_struct(imag_part);
    init_mpfr_value(imag_part, prec);
    mpfr_set(imag_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
  }
}


inline mp_prec_t cxfloat::get_prec() const
{
  return mpfr_get_prec(real_part);
}


inline void cxfloat::get_real(refloat& val) const
{
  mpfr_set(val.value, real_part, GMP_RNDN);
}


inline void cxfloat::get_imag(refloat& val) const 
{
  mpfr_set(val.value, imag_part, GMP_RNDN);
}


inline void swap(cxfloat& a, cxfloat& b)
{
  std::swap( a.real_part, b.real_part );
  std::swap( a.imag_part, b.imag_part );
}


inline void cxfloat::neg()
{
  if(real_part != nullptr){mpfr_neg(real_part, real_part, GMP_RNDN);}
  if(imag_part != nullptr){mpfr_neg(imag_part, imag_part, GMP_RNDN);}
}


inline void neg(const cxfloat& a, cxfloat& c)
{
  c = a;
  c.neg();
}


inline void cxfloat::add(const cxfloat& x)
{
  if( (real_part != nullptr) && (x.real_part != nullptr) ){mpfr_add(real_part, real_part, x.real_part, GMP_RNDN);}
  else if( (real_part == nullptr) && (x.real_part != nullptr) ){set_real(x);}
  if( (imag_part != nullptr) && (x.imag_part != nullptr) ){mpfr_add(imag_part, imag_part, x.imag_part, GMP_RNDN);}
  else if( (imag_part == nullptr) && (x.imag_part != nullptr) ){set_imag(x);}
}

inline void cxfloat::add(const refloat& x)
{
  if( (real_part != nullptr) && (x.value != nullptr) ){mpfr_add(real_part, real_part, x.value, GMP_RNDN);}
}


inline void add(const cxfloat& a, const cxfloat& b, cxfloat& c)
{
  c = a;
  c.add(b);
}


inline void cxfloat::sub(const cxfloat& x)
{
  if( (real_part != nullptr) && (x.real_part != nullptr) ){mpfr_sub(real_part, real_part, x.real_part, GMP_RNDN);}
  else if( (real_part == nullptr) && (x.real_part != nullptr) ){set_real(x);mpfr_neg(real_part, real_part, GMP_RNDN);}
  if( (imag_part != nullptr) && (x.imag_part != nullptr) ){mpfr_sub(imag_part, imag_part, x.imag_part, GMP_RNDN);}
  else if( (imag_part == nullptr) && (x.imag_part != nullptr) ){set_imag(x);mpfr_neg(imag_part, imag_part, GMP_RNDN);}
}

inline void cxfloat::sub(const refloat& x)
{
  if( (real_part != nullptr) && (x.value != nullptr) ){mpfr_sub(real_part, real_part, x.value, GMP_RNDN);}
}


inline void sub(const cxfloat& a, const cxfloat& b, cxfloat& c)
{
  c = a;
  c.sub(b);
}


inline void cxfloat::mul(const cxfloat& x)
{
  if( (real_part != nullptr) )
  {
    if( imag_part != nullptr )
    {
      if( x.real_part != nullptr )
      {
        if(x.imag_part != nullptr)
        {
          if(mpfr_get_prec(x.real_part) >= (1<<11))
          {
            mpfr_add((*cxfloat_cache1).real_part, real_part, imag_part, GMP_RNDN);
            mpfr_sub((*cxfloat_cache1).imag_part, x.imag_part, x.real_part, GMP_RNDN);
            mpfr_mul((*cxfloat_cache1).real_part, (*cxfloat_cache1).real_part, x.real_part, GMP_RNDN);
            mpfr_mul((*cxfloat_cache1).imag_part, (*cxfloat_cache1).imag_part, real_part, GMP_RNDN);
            mpfr_add(real_part, x.real_part, x.imag_part, GMP_RNDN);
            mpfr_mul(real_part, real_part, imag_part, GMP_RNDN);
            mpfr_sub(real_part, (*cxfloat_cache1).real_part, real_part, GMP_RNDN);
            mpfr_add(imag_part, (*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
          }
          else
          {
            mpfr_mul((*cxfloat_cache1).real_part, real_part, x.real_part, GMP_RNDN);
            mpfr_mul((*cxfloat_cache1).imag_part, real_part, x.imag_part, GMP_RNDN);
            mpfr_mul(real_part, imag_part, x.imag_part, GMP_RNDN);
            mpfr_neg(real_part, real_part, GMP_RNDN);
            mpfr_add(real_part, real_part, (*cxfloat_cache1).real_part, GMP_RNDN);
            mpfr_mul(imag_part, imag_part, x.real_part, GMP_RNDN);
            mpfr_add(imag_part, imag_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
          }
        }
        else
        {
          mpfr_mul(real_part, real_part, x.real_part, GMP_RNDN);
          mpfr_mul(imag_part, imag_part, x.real_part, GMP_RNDN);
        }
      }
      else
      {
        if(x.imag_part != nullptr)
        {
          mpfr_mul((*cxfloat_cache1).real_part, imag_part, x.imag_part, GMP_RNDN);
          mpfr_mul(imag_part, real_part, x.imag_part, GMP_RNDN);
          mpfr_neg(real_part, (*cxfloat_cache1).real_part, GMP_RNDN);
        }
        else
        {
          clear();
        }
      }
    }
    else
    {
      if(x.real_part != nullptr)
      {
        if(x.imag_part != nullptr)
        {
          set_imag(0);
          mpfr_mul(imag_part, real_part, x.imag_part, GMP_RNDN);
          mpfr_mul(real_part, real_part, x.real_part, GMP_RNDN);
        }
        else
        {
          mpfr_mul(real_part, real_part, x.real_part, GMP_RNDN);
        }
      }
      else
      {
        if(x.imag_part != nullptr)
        {
          mpfr_mul(real_part, real_part, x.imag_part, GMP_RNDN);
          std::swap(real_part, imag_part);
        }
        else
        {
          clear();
        }
      }
    }
  }
  else
  {
    if( imag_part != nullptr )
    {
      if( x.real_part != nullptr )
      {
        if(x.imag_part != nullptr)
        {
          set_real(0);
          mpfr_mul(real_part, imag_part, x.imag_part, GMP_RNDN);
          mpfr_neg(real_part, real_part, GMP_RNDN);
          mpfr_mul(imag_part, imag_part, x.real_part, GMP_RNDN);
        }
        else
        {
          mpfr_mul(imag_part, imag_part, x.real_part, GMP_RNDN);
        }
      }
      else
      {
        if(x.imag_part != nullptr)
        {
          mpfr_mul(imag_part, imag_part, x.imag_part, GMP_RNDN);
          std::swap(real_part, imag_part);
          mpfr_neg(real_part, real_part, GMP_RNDN);
        }
        else
        {
          clear();
        }
      }
    }
    else
    {
      clear();
    }
  }
}


inline void cxfloat::mul(const refloat& x)
{
  mpfr_mul(real_part, real_part, x.value, GMP_RNDN);
  mpfr_mul(imag_part, imag_part, x.value, GMP_RNDN);
}


inline void mul(const cxfloat& a, const cxfloat& b, cxfloat& c)
{
  c = a;
  c.mul(b);
}


inline void mul(const refloat& a, const cxfloat& b, cxfloat& c)
{
  c = b;
  c.mul(a);
}


inline void cxfloat::sqr()
{
  if( (real_part != nullptr) && (real_part != nullptr) )
  {
    mpfr_sqr((*cxfloat_cache1).real_part, real_part, GMP_RNDN);
    mpfr_sqr((*cxfloat_cache1).imag_part, imag_part, GMP_RNDN);
    mpfr_neg((*cxfloat_cache1).imag_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
    mpfr_add((*cxfloat_cache1).real_part, (*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
    mpfr_mul((*cxfloat_cache1).imag_part, real_part, imag_part, GMP_RNDN);
    mpfr_mul_2ui(imag_part, (*cxfloat_cache1).imag_part, 1, GMP_RNDN);
    mpfr_set(real_part, (*cxfloat_cache1).real_part, GMP_RNDN);
  }
  else if( (real_part != nullptr) )
  {
    mpfr_sqr(real_part, real_part, GMP_RNDN);
  }
  else if( (imag_part != nullptr) )
  {
    mpfr_sqr(real_part, imag_part, GMP_RNDN);
    mpfr_neg(real_part, real_part, GMP_RNDN);
  }
}


inline void sqr(const cxfloat& a, cxfloat& c)
{
  c = a;
  c.sqr();
}


inline void cxfloat::div(const cxfloat& x)
{
  (*cxfloat_cache2).set(x);
  mpfr_neg((*cxfloat_cache2).imag_part, (*cxfloat_cache2).imag_part, GMP_RNDN);
  this->mul((*cxfloat_cache2));
  
  if( (x.real_part != nullptr) && (x.real_part != nullptr) )
  {
    mpfr_sqr((*cxfloat_cache1).real_part, x.real_part, GMP_RNDN);
    mpfr_sqr((*cxfloat_cache1).imag_part, x.imag_part, GMP_RNDN);
    mpfr_add((*cxfloat_cache1).real_part, (*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
  }
  else if(x.real_part != nullptr)
  {
    mpfr_sqr((*cxfloat_cache1).real_part, x.real_part, GMP_RNDN);
  }
  else if(x.imag_part != nullptr)
  {
    mpfr_sqr((*cxfloat_cache1).real_part, x.imag_part, GMP_RNDN);
  }
  
  if(real_part != nullptr){mpfr_div(real_part, real_part, (*cxfloat_cache1).real_part, GMP_RNDN);}
  if(imag_part != nullptr){mpfr_div(imag_part, imag_part, (*cxfloat_cache1).real_part, GMP_RNDN);}
}


inline void cxfloat::div(const refloat& x)
{
  if(real_part != nullptr){mpfr_div(real_part, real_part, x.value, GMP_RNDN);}
  if(imag_part!= nullptr){mpfr_div(imag_part, imag_part, x.value, GMP_RNDN);}
}


inline void div(const cxfloat& a, const cxfloat& b, cxfloat& c)
{
  c.set(a);
  c.div(b);
}


inline void cxfloat::sqrt()
{
  if(real_part != nullptr)
  {
    if(imag_part != nullptr)
    {
      if(mpfr_zero_p(imag_part)!=0)
      {
        mpfr_sqrt(real_part, real_part, GMP_RNDN);
      }
      else
      {
        mpfr_sqr((*cxfloat_cache1).real_part, real_part, GMP_RNDN);
        mpfr_sqr((*cxfloat_cache1).imag_part, imag_part, GMP_RNDN);
        mpfr_add((*cxfloat_cache1).real_part, (*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
        mpfr_sqrt((*cxfloat_cache1).real_part, (*cxfloat_cache1).real_part, GMP_RNDN);
        mpfr_add(real_part, real_part, (*cxfloat_cache1).real_part, GMP_RNDN);
        mpfr_div_2si(real_part, real_part, 1, GMP_RNDN);
        mpfr_sqrt(real_part, real_part, GMP_RNDN);
        mpfr_div(imag_part, imag_part, real_part, GMP_RNDN);
        mpfr_div_2si(imag_part, imag_part, 1, GMP_RNDN);
      }
    }
    else
    {
      if( mpfr_cmp_si(real_part, 0) < 0 )
      {
        std::swap(real_part, imag_part);
        mpfr_neg(imag_part, imag_part, GMP_RNDN);
        mpfr_sqrt(imag_part, imag_part, GMP_RNDN);
      }
      else
      {
        mpfr_sqrt((*cxfloat_cache1).real_part, real_part, GMP_RNDN);
      }
    }
  }
  else
  {
    if(imag_part != nullptr)
    {
      set_real(0);
      sqrt();
    }
  }
}


inline void sqrt(const cxfloat& a, cxfloat& b)
{
  b.set(a);
  b.sqrt();
}


inline void cxfloat::mul_power_2(long int p)
{
  if(p>0)
  {
    if(real_part != nullptr){mpfr_mul_2si(real_part, real_part, p, GMP_RNDN);}
    if(imag_part != nullptr){mpfr_mul_2si(imag_part, imag_part, p, GMP_RNDN);}
  }
  else
  {
    if(real_part != nullptr){mpfr_div_2si(real_part, real_part, -p, GMP_RNDN);}
    if(imag_part!= nullptr){mpfr_div_2si(imag_part, imag_part, -p, GMP_RNDN);}
  }
}


inline void mul_power_2(const cxfloat& a, long int p, cxfloat& c)
{
  c.set(a);
  c.mul_power_2(p);
}


inline void cxfloat::mul_by_i()
{
  std::swap(real_part, imag_part);
  if(real_part != nullptr){mpfr_neg(real_part, real_part, GMP_RNDN);}
}


inline void mul_by_i(const cxfloat& a, cxfloat& b)
{
  b.set(a);
  b.mul_by_i();
}


inline void cxfloat::conj()
{
  if(imag_part!= nullptr){mpfr_neg(imag_part, imag_part, GMP_RNDN);}
}


inline void conj(const cxfloat& a, cxfloat& b)
{
  b.set(a);
  b.conj();
}


inline void cxfloat::reciprocal()
{
  if( (real_part != nullptr) && (imag_part != nullptr) )
  {
    mpfr_sqr((*cxfloat_cache1).real_part, real_part, GMP_RNDN);
    mpfr_sqr((*cxfloat_cache1).imag_part, imag_part, GMP_RNDN);
    mpfr_add((*cxfloat_cache1).real_part, (*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
    mpfr_neg(imag_part, imag_part, GMP_RNDN);
    mpfr_div(real_part, real_part, (*cxfloat_cache1).real_part, GMP_RNDN);
    mpfr_div(imag_part, imag_part, (*cxfloat_cache1).real_part, GMP_RNDN);
  }
  else if( real_part != nullptr )
  {
    mpfr_set_si((*cxfloat_cache1).real_part, 1, GMP_RNDN);
    mpfr_div(real_part, (*cxfloat_cache1).real_part, real_part, GMP_RNDN);
  }
  else if( imag_part != nullptr )
  {
    mpfr_set_si((*cxfloat_cache1).imag_part, 1, GMP_RNDN);
    mpfr_div(imag_part, (*cxfloat_cache1).imag_part, imag_part, GMP_RNDN);
    mpfr_neg(imag_part, imag_part, GMP_RNDN);
  }
}


inline void magnitude(const cxfloat& a, cxfloat& b)
{
  b.clear_imag();
  if(b.real_part == nullptr){b.set_real(0);}
  if( (a.real_part != nullptr) && (a.imag_part != nullptr) )
  {
    mpfr_sqr((*cxfloat_cache1).real_part, a.real_part, GMP_RNDN);
    mpfr_sqr((*cxfloat_cache1).imag_part, a.imag_part, GMP_RNDN);
    mpfr_add((*cxfloat_cache1).real_part, (*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
    mpfr_sqrt(b.real_part, (*cxfloat_cache1).real_part, GMP_RNDN);
  }
  else if( a.real_part != nullptr )
  {
    mpfr_set(b.real_part, a.real_part, GMP_RNDN);
  }
  else if( a.imag_part != nullptr )
  {
    mpfr_set(b.real_part, a.imag_part, GMP_RNDN);
  }
  else
  {
    mpfr_set_si(b.real_part, 0, GMP_RNDN);
  }
}


inline void magnitude(const cxfloat& a, refloat& b)
{
  if( (a.real_part != nullptr) && (a.imag_part != nullptr) )
  {
    mpfr_sqr((*cxfloat_cache1).real_part, a.real_part, GMP_RNDN);
    mpfr_sqr((*cxfloat_cache1).imag_part, a.imag_part, GMP_RNDN);
    mpfr_add((*cxfloat_cache1).real_part, (*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
    mpfr_sqrt(b.value, (*cxfloat_cache1).real_part, GMP_RNDN);
  }
  else if( a.real_part != nullptr )
  {
    mpfr_set(b.value, a.real_part, GMP_RNDN);
  }
  else if( a.imag_part != nullptr )
  {
    mpfr_set(b.value, a.imag_part, GMP_RNDN);
  }
  else
  {
    mpfr_set_si(b.value, 0, GMP_RNDN);
  }
}


inline void sin_cos(const cxfloat& z, cxfloat& sn, cxfloat& cs)
{
  if( (z.real_part != nullptr) && (z.imag_part != nullptr) )
  {
    if( sn.real_part == nullptr ){sn.set_real(0);}
    if( cs.real_part == nullptr ){cs.set_real(0);}
    if( sn.imag_part == nullptr ){sn.set_imag(0);}
    if( cs.imag_part == nullptr ){cs.set_imag(0);}
    
    mpfr_sin_cos((*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, z.real_part, GMP_RNDN);
    if(mpfr_zero_p(z.imag_part)!=0)
    {
      mpfr_set(sn.real_part, (*cxfloat_cache1).real_part, GMP_RNDN);
      mpfr_set(cs.real_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
      mpfr_set_si(sn.imag_part, 0, GMP_RNDN);
      mpfr_set_si(cs.imag_part, 0, GMP_RNDN);
    }
    else
    {
      mpfr_sinh_cosh((*cxfloat_cache2).real_part, (*cxfloat_cache2).imag_part, z.imag_part, GMP_RNDN);
      mpfr_mul(sn.real_part, (*cxfloat_cache1).real_part, (*cxfloat_cache2).imag_part, GMP_RNDN);
      mpfr_mul(sn.imag_part, (*cxfloat_cache1).imag_part, (*cxfloat_cache2).real_part, GMP_RNDN);
      mpfr_mul(cs.real_part, (*cxfloat_cache1).imag_part, (*cxfloat_cache2).imag_part, GMP_RNDN);
      mpfr_mul(cs.imag_part, (*cxfloat_cache1).real_part, (*cxfloat_cache2).real_part, GMP_RNDN);
      mpfr_neg(cs.imag_part, cs.imag_part, GMP_RNDN);
    }
  }
  else if( z.real_part != nullptr )
  {
    if( sn.real_part == nullptr ){sn.set_real(0);}
    if( cs.real_part == nullptr ){cs.set_real(0);}
    mpfr_sin_cos((*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, z.real_part, GMP_RNDN);
    mpfr_set(sn.real_part, (*cxfloat_cache1).real_part, GMP_RNDN);
    mpfr_set(cs.real_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
    sn.clear_imag();
    cs.clear_imag();
  }
  else if( z.imag_part != nullptr )
  {
    if( sn.imag_part == nullptr ){sn.set_imag(0);}
    if( cs.real_part == nullptr ){cs.set_real(0);}
    mpfr_sinh_cosh((*cxfloat_cache2).real_part, (*cxfloat_cache2).imag_part, z.imag_part, GMP_RNDN);
    sn.clear_real();
    mpfr_set(sn.imag_part, (*cxfloat_cache2).real_part, GMP_RNDN);
    mpfr_set(cs.real_part, (*cxfloat_cache2).imag_part, GMP_RNDN);
    cs.clear_imag();
  }
}


inline void cos(const cxfloat& z, cxfloat& cs)
{
  if( (z.real_part != nullptr) && (z.imag_part != nullptr) )
  {
    if( cs.real_part == nullptr ){cs.set_real(0);}
    if( cs.imag_part == nullptr ){cs.set_imag(0);}
    
    if(mpfr_zero_p(z.imag_part)!=0)
    {
      mpfr_cos(cs.real_part, z.real_part, GMP_RNDN);
      mpfr_set_si(cs.imag_part, 0, GMP_RNDN);
    }
    else
    {
      mpfr_sin_cos((*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, z.real_part, GMP_RNDN);
      mpfr_sinh_cosh((*cxfloat_cache2).real_part, (*cxfloat_cache2).imag_part, z.imag_part, GMP_RNDN);
      mpfr_mul(cs.real_part, (*cxfloat_cache1).imag_part, (*cxfloat_cache2).imag_part, GMP_RNDN);
      mpfr_mul(cs.imag_part, (*cxfloat_cache1).real_part, (*cxfloat_cache2).real_part, GMP_RNDN);
      mpfr_neg(cs.imag_part, cs.imag_part, GMP_RNDN);
    }
  }
  else if( z.real_part != nullptr )
  {
    if( cs.real_part == nullptr ){cs.set_real(0);}
    mpfr_sin_cos((*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, z.real_part, GMP_RNDN);
    mpfr_set(cs.real_part, (*cxfloat_cache1).imag_part, GMP_RNDN);
    cs.clear_imag();
  }
  else if( z.imag_part != nullptr )
  {
    if( cs.real_part == nullptr ){cs.set_real(0);}
    mpfr_sinh_cosh((*cxfloat_cache2).real_part, (*cxfloat_cache2).imag_part, z.imag_part, GMP_RNDN);
    mpfr_set(cs.real_part, (*cxfloat_cache2).imag_part, GMP_RNDN);
    cs.clear_imag();
  }
}


inline void sin(const cxfloat& z, cxfloat& sn)
{
  if( (z.real_part != nullptr) && (z.imag_part != nullptr) )
  {
    if( sn.real_part == nullptr ){sn.set_real(0);}
    if( sn.imag_part == nullptr ){sn.set_imag(0);}
    if(mpfr_zero_p(z.imag_part)!=0)
    {
      mpfr_sin(sn.real_part, z.real_part, GMP_RNDN);
      mpfr_set_si(sn.imag_part, 0, GMP_RNDN);
    }
    else
    {
      mpfr_sin_cos((*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, z.real_part, GMP_RNDN);
      mpfr_sinh_cosh((*cxfloat_cache2).real_part, (*cxfloat_cache2).imag_part, z.imag_part, GMP_RNDN);
      mpfr_mul(sn.real_part, (*cxfloat_cache1).real_part, (*cxfloat_cache2).imag_part, GMP_RNDN);
      mpfr_mul(sn.imag_part, (*cxfloat_cache1).imag_part, (*cxfloat_cache2).real_part, GMP_RNDN);
    }
  }
  else if( z.real_part != nullptr )
  {
    if( sn.real_part == nullptr ){sn.set_real(0);}
    mpfr_sin_cos((*cxfloat_cache1).real_part, (*cxfloat_cache1).imag_part, z.real_part, GMP_RNDN);
    mpfr_set(sn.real_part, (*cxfloat_cache1).real_part, GMP_RNDN);
    sn.clear_imag();
  }
  else if( z.imag_part != nullptr )
  {
    if( sn.imag_part == nullptr ){sn.set_imag(0);}
    mpfr_sinh_cosh((*cxfloat_cache2).real_part, (*cxfloat_cache2).imag_part, z.imag_part, GMP_RNDN);
    sn.clear_real();
    mpfr_set(sn.imag_part, (*cxfloat_cache2).real_part, GMP_RNDN);
  }
}


inline void exp(const cxfloat& z, cxfloat& e)
{
  if( (z.real_part != nullptr) && (z.imag_part != nullptr) )
  {
    if(e.real_part == nullptr){e.set_real(0);}
    if(e.imag_part == nullptr){e.set_imag(0);}
    
    int rez = mpfr_zero_p(z.real_part);
    int imz = mpfr_zero_p(z.imag_part);
    
    if(imz==0)
    {
      mpfr_sin_cos(e.imag_part, e.real_part, z.imag_part, GMP_RNDN);
    }
    else
    {
      e.set(1,0);
    }
    if(rez==0)
    {
      mpfr_exp((*cxfloat_cache1).real_part, z.real_part, GMP_RNDN);
      mpfr_mul(e.real_part, (*cxfloat_cache1).real_part, e.real_part, GMP_RNDN);
      mpfr_mul(e.imag_part, (*cxfloat_cache1).real_part, e.imag_part, GMP_RNDN);
    }
  }
  else if(z.real_part != nullptr)
  {
    if(e.real_part == nullptr){e.set_real(0);}
    e.clear_imag();
    mpfr_exp(e.real_part, z.real_part, GMP_RNDN);
  }
  else if(z.imag_part != nullptr)
  {
    if(e.real_part == nullptr){e.set_real(0);}
    if(e.imag_part == nullptr){e.set_imag(0);}
    mpfr_sin_cos(e.imag_part, e.real_part, z.imag_part, GMP_RNDN);
  }
}


inline void log(const cxfloat& z, cxfloat& l)
{
  if( (z.real_part != nullptr) && (z.imag_part != nullptr) )
  {
    if(l.real_part == nullptr){l.set_real(0);}
    if(l.imag_part == nullptr){l.set_imag(0);}
    
    magnitude( z, *cxfloat_cache2 );
    mpfr_log( l.real_part, cxfloat_cache2->real_part, GMP_RNDN );
    mpfr_atan2( l.imag_part, z.imag_part, z.real_part, GMP_RNDN );
    cxfloat_cache2->set_imag(0);
  }
  else if(z.real_part != nullptr)
  {
    if(l.real_part == nullptr){l.set_real(0);}
    l.clear_imag();
    mpfr_log( l.real_part, z.real_part, GMP_RNDN );
  }
  else if(z.imag_part != nullptr)
  {
    if(l.real_part == nullptr){l.set_real(0);}
    if(l.imag_part == nullptr){l.set_imag(0);}
    mpfr_log( l.real_part, z.imag_part, GMP_RNDN );
    mpfr_const_pi(l.imag_part, GMP_RNDN);
    mpfr_div_2si(l.imag_part, l.imag_part, 1, GMP_RNDN);
    if( mpfr_cmp_si(z.imag_part, 0) < 0 ){ mpfr_neg(l.imag_part, l.imag_part, GMP_RNDN); }
  }
}


inline void pow(unsigned long int n, const cxfloat& x, cxfloat& v)
{
  if(x.real_part != nullptr){v.set(1, 0, mpfr_get_prec(x.real_part));}
  if(x.imag_part != nullptr){v.set(1, 0, mpfr_get_prec(x.imag_part));}
  else{return;}
  (*cxfloat_cache2).set(x);
  
  while(n>0)
  {
    if(n&1)
    {
      v.mul((*cxfloat_cache2));
    }
    (*cxfloat_cache2).sqr();
    n = (n>>1);
  }
  if(cxfloat_cache2->real_part == nullptr){cxfloat_cache2->set_real(0);}
  if(cxfloat_cache2->imag_part == nullptr){cxfloat_cache2->set_imag(0);}
}


inline int compare_real(const cxfloat& a, const cxfloat& b)
{
  return mpfr_cmp(a.real_part, b.real_part);
}


inline int compare_real_abs(const cxfloat& a, const cxfloat& b)
{
  return mpfr_cmpabs(a.real_part, b.real_part);
}


inline int compare_imag(const cxfloat& a, const cxfloat& b)
{
  return mpfr_cmp(a.imag_part, b.imag_part);
}


inline int compare_imag_abs(const cxfloat& a, const cxfloat& b)
{
  return mpfr_cmpabs(a.imag_part, b.imag_part);
}









#endif
